MODULE dcmip_baro_wave_test_mod
  
  use const_mod, only: r8
  use formula_mod
  use latlon_parallel_mod
  use block_mod
  use operators_mod
  !=======================================================================
  !
  !  Date:  July 29, 2015
  !
  !  Functions for setting up idealized initial conditions for the
  !  Ullrich, Melvin, Staniforth and Jablonowski baroclinic instability.
  !
  !  SUBROUTINE baroclinic_wave_sample(
  !    deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)
  !
  !  Options:
  !     deep    deep atmosphere (1 = yes or 0 = no)
  !    moist    include moisture (1 = yes or 0 = no)
  !    pertt    type of perturbation (0 = exponential, 1 = stream function)
  !        X    Earth scaling factor
  !
  !  Given a point specified by: 
  !      lon    longitude (radians) 
  !      lat    latitude (radians) 
  !      p/z    pressure (Pa) / height (m)
  !  zcoords    1 if z is specified, 0 if p is specified
  !
  !  the functions will return:
  !        p    pressure if z is specified and zcoords = 1 (Pa)
  !        u    zonal wind (m s^-1)
  !        v    meridional wind (m s^-1)
  !        t    temperature (K)
  !   thetav    virtual potential temperature (K)
  !     phis    surface geopotential (m^2 s^-2)
  !       ps    surface pressure (Pa)
  !      rho    density (kj m^-3)
  !        q    water vapor mixing ratio (kg/kg)
  !
  !
  !  Author: Paul Ullrich
  !          University of California, Davis
  !          Email: paullrich@ucdavis.edu
  !
  !=======================================================================
  
    IMPLICIT NONE

    private
  
    public deepatm_baro_wave_test_set_ic
  
  !=======================================================================
  !    Physical constants
  !=======================================================================
  
    REAL(8), PARAMETER ::               &
         a     = 6371220.0d0,           & ! Reference Earth's Radius (m)
         Rd    = 287.0d0,               & ! Ideal gas const dry air (J kg^-1 K^1)
         g     = 9.80616d0,             & ! Gravity (m s^2)
         cp    = 1004.5d0,              & ! Specific heat capacity (J kg^-1 K^1)
         Lvap  = 2.5d6,                 & ! Latent heat of vaporization of water
         Rvap  = 461.5d0,               & ! Ideal gas constnat for water vapor
         Mvap  = 0.608d0,               & ! Ratio of molar mass of dry air/water
         pi    = 3.14159265358979d0,    & ! pi
         p0    = 100000.0d0,            & ! surface pressure (Pa)
         kappa = 2.d0/7.d0,             & ! Ratio of Rd to cp
         omega = 7.29212d-5,            & ! Reference rotation rate of the Earth (s^-1)
         deg2rad  = pi/180.d0             ! Conversion factor of degrees to radians
  
  !=======================================================================
  !    Test case parameters
  !=======================================================================
    REAL(8), PARAMETER ::               &
         T0E        = 310.d0     ,      & ! temperature at equatorial surface (K)
         T0P        = 240.d0     ,      & ! temperature at polar surface (K)
         B          = 2.d0       ,      & ! jet half-width parameter
         K          = 3.d0       ,      & ! jet width parameter
         lapse      = 0.005d0             ! lapse rate parameter
  
    REAL(8), PARAMETER ::               &
         pertu0     = 0.5d0      ,      & ! SF Perturbation wind velocity (m/s)
         pertr      = 1.d0/6.d0  ,      & ! SF Perturbation radius (Earth radii)
         pertup     = 1.0d0      ,      & ! Exp. perturbation wind velocity (m/s)
         pertexpr   = 0.1d0      ,      & ! Exp. perturbation radius (Earth radii)
         pertlon    = pi/9.d0    ,      & ! Perturbation longitude
         pertlat    = 2.d0*pi/9.d0,     & ! Perturbation latitude
         pertz      = 15000.d0   ,      & ! Perturbation height cap
         dxepsilon  = 1.d-5               ! Small value for numerical derivatives
   
    REAL(8), PARAMETER ::               &
         moistqlat  = 2.d0*pi/9.d0,     & ! Humidity latitudinal width
         moistqp    = 34000.d0,         & ! Humidity vertical pressure width
         moisttr    = 0.1d0,            & ! Vertical cut-off pressure for humidity
         moistqs    = 1.d-12,           & ! Humidity above cut-off
         moistq0    = 0.018d0,          & ! Maximum specific humidity
         moistqr    = 0.9d0,            & ! Maximum saturation ratio
         moisteps   = 0.622d0,          & ! Ratio of gas constants
         moistT0    = 273.16d0,         & ! Reference temperature (K)
         moistE0Ast = 610.78d0            ! Saturation vapor pressure at T0 (Pa) 
  
  CONTAINS


  subroutine deepatm_baro_wave_test_set_ic(block)
  
    type(block_type), intent(inout), target :: block

    ! real(r8) etav, eta, tbar, ubar, gzbar, sin_lat, cos_lat, lamda, half_lon, d, phi, big_u
    integer :: deep,moist,pertt,zcoords

    integer i, j, k

    real(r8) :: X
    real(r8) :: gm_lon, gm_lat
    real(r8) :: rho_gm, q_gm, u_gm, v_gm, gzs_gm, t_gm, pt_gm, ps_gm ! not used variable
    
    ! real(r8) :: rp !! calc r from p
    ! real(r8),allocatable :: r_full_3d, r_half_3d
    associate (mesh   => block%mesh            , &
               u      => block%dstate(1)%u_lon , &
               v      => block%dstate(1)%v_lat , &
               mgs    => block%dstate(1)%mgs   , &
               mg     => block%dstate(1)%mg    , &
               mg_lev => block%dstate(1)%mg_lev, &
               t      => block%aux%t           , &
               pt     => block%dstate(1)%pt    , &
               !rho    => block%dstate(1)%rhod  , &
               gz_lev => block%dstate(1)%gz_lev, &
               gz     => block%dstate(1)%gz    , &
               gzs    => block%static%gzs)
    mgs%d = 1.0e5_r8
    ! v  %d = 0

    deep    = 1
    moist   = 0
    pertt   = 1
    zcoords = 0 
    X       = 1.0_r8
    ! allocate(r_full_3d(mesh%full_ids:mesh%full_ide,mesh%full_ids:mesh%full_ids,mesh%full_ids:mesh%full_ids))
    call calc_mg(block, block%dstate(1))

    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde
        gm_lat = mesh%full_lat(j)
        do i = mesh%full_ids, mesh%full_ide
          gm_lon = mesh%full_lon(i) 
          ! determine the temperature, height
          call baroclinic_wave_test(deep    = deep        , &
            &                       moist   = moist       , &
            &                       pertt   = pertt       , &
            &                       X       = X           , &
            &                       lon     = gm_lon      , &
            &                       lat     = gm_lat      , &
            &                       p       = mg%d(i,j,k) , &
            &                       z       = gz%d(i,j,k) , & 
            &                       zcoords = zcoords     , &
            &                       u       = u_gm        , &  ! not use , because it is at half_id
            &                       v       = v_gm        , &  ! not use , because it is at half_jd
            &                       t       = t%d(i,j,k)  , &  ! out
            &                       thetav  = pt%d(i,j,k) , &  ! out
            &                       phis    = gzs_gm      , &  ! out
            &                       ps      = ps_gm       , &
            &                       rho     = rho_gm      , &
            &                       q       = q_gm        )
            pt%d(i,j,k) = modified_potential_temperature(t%d(i,j,k), mg%d(i,j,k), 0.0_r8)
        end do
      end do
    end do
    call fill_halo(pt)

    do k = mesh%half_kds, mesh%half_kde
      do j = mesh%full_jds, mesh%full_jde
        gm_lat = mesh%full_lat(j)
        do i = mesh%full_ids, mesh%full_ide
          gm_lon = mesh%full_lon(i) 
          call baroclinic_wave_test(deep    = deep            , &
            &                       moist   = moist           , &
            &                       pertt   = pertt           , &
            &                       X       = X               , &
            &                       lon     = gm_lon          , &
            &                       lat     = gm_lat          , &
            &                       p       = mg_lev%d(i,j,k) , &
            &                       z       = gz_lev%d(i,j,k) , & 
            &                       zcoords = zcoords         , &
            &                       u       = u_gm            , &  ! not use , because it is at half_id
            &                       v       = v_gm            , &  ! not use , because it is at half_jd
            &                       t       = t_gm            , &  ! out
            &                       thetav  = pt_gm           , &  ! out
            &                       phis    = gzs_gm          , &  ! out
            &                       ps      = ps_gm           , &
            &                       rho     = rho_gm          , &
            &                       q       = q_gm            )
        end do
      end do
    end do

    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde
        gm_lat = mesh%full_lat(j)
        do i = mesh%half_ids, mesh%half_ide
          gm_lon = mesh%half_lon(i) 

          call init_velocity_field(deep = deep               , & !in
          &                        pertt= pertt              , & !in
          &                        X    = X                  , & !in
          &                        lon  = gm_lon             , & !in
          &                        lat  = gm_lat             , & !in
          &                        z    = gz%d(i,j,k)        , & !in
          &                        tempv= t%d(i,j,k)         , & !in
          &                        u    = u%d(i,j,k)         , & !out
          &                        v    = v_gm                 ) !out)
        end do
      end do

      do j = mesh%half_jds, mesh%half_jde
        gm_lat = mesh%half_lat(j)
        do i = mesh%full_ids, mesh%full_ide
          gm_lon = mesh%full_lon(i) 

          call init_velocity_field(deep = deep               , & !in
          &                        pertt= pertt              , & !in
          &                        X    = X                  , & !in
          &                        lon  = gm_lon             , & !in
          &                        lat  = gm_lat             , & !in
          &                        z    = gz%d(i,j,k)        , & !in
          &                        tempv= t%d(i,j,k)         , & !in
          &                        u    = u_gm               , & !out
          &                        v    = v%d(i,j,k)         ) !out)
        end do
      end do
    end do
    call fill_halo(u)
    call fill_halo(v)

    gz%d     = gz%d     * g
    gz_lev%d = gz_lev%d * g
    gzs%d    = 0

    call fill_halo(gz)
    call fill_halo(gz_lev)
    call fill_halo(gzs)
    
    end associate

  end subroutine deepatm_baro_wave_test_set_ic
  
  !=======================================================================
  !    Generate the baroclinic instability initial conditions   
  !=======================================================================
  ! ICON has a new way to realise the baro_deep/shallow wave, i follow that
    SUBROUTINE baroclinic_wave_test(deep,moist,pertt,X,lon,lat,p,z,zcoords,u,v,t,thetav,phis,ps,rho,q) !&
      !BIND(c, name = "baroclinic_wave_test")
   
      !IMPLICIT NONE
  
  !-----------------------------------------------------------------------
  !     input/output params parameters at given location
  !-----------------------------------------------------------------------
      INTEGER, INTENT(IN)  :: &
                  deep,       & ! Deep (1) or Shallow (0) test case
                  moist,      & ! Moist (1) or Dry (0) test case
                  pertt         ! Perturbation type
  
      REAL(r8), INTENT(IN)  :: &
                  lon,        & ! Longitude (radians)
                  lat,        & ! Latitude (radians)
                  X             ! Earth scaling parameter
  
      REAL(r8), INTENT(INOUT) :: &
                  p,            & ! Pressure (Pa)
                  z               ! Altitude (m)
  
      INTEGER, INTENT(IN) :: zcoords     ! 1 if z coordinates are specified
                                         ! 0 if p coordinates are specified
  
      REAL(r8), INTENT(OUT) :: &
                  u,          & ! Zonal wind (m s^-1)
                  v,          & ! Meridional wind (m s^-1)
                  t,          & ! Temperature (K)
                  thetav,     & ! Virtual potential temperature (K)
                  phis,       & ! Surface Geopotential (m^2 s^-2)
                  ps,         & ! Surface Pressure (Pa)
                  rho,        & ! density (kg m^-3)
                  q             ! water vapor mixing ratio (kg/kg)
  
      !------------------------------------------------
      !   Local variables
      !------------------------------------------------
      ! REAL(8) :: aref, omegaref
      ! REAL(8) :: T0, constH, constC, scaledZ, inttau2, rratio
      ! REAL(8) :: inttermU, bigU, rcoslat, omegarcoslat
      ! REAL(8) :: eta, qratio, qnum, qden
      REAL(r8) :: eta
      !------------------------------------------------
      !   Pressure and temperature
      !------------------------------------------------
      if (zcoords .eq. 1) then
        CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
      else
        CALL evaluate_z_temperature(deep, X, lon, lat, p, z, t)
      end if
  
      ! !------------------------------------------------
      ! !   Compute test case constants
      ! !------------------------------------------------
      ! aref = a / X
      ! omegaref = omega * X
  
      ! T0 = 0.5d0 * (T0E + T0P)
  
      ! constH = Rd * T0 / g
  
      ! constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
  
      ! scaledZ = z / (B * constH)
  
      ! inttau2 = constC * z * exp(- scaledZ**2)
  
      ! ! radius ratio
      ! if (deep .eq. 0) then
      !   rratio = 1.d0
      ! else
      !   rratio = (z + aref) / aref;
      ! end if
  
      ! !-----------------------------------------------------
      ! !   Initialize surface pressure
      ! !-----------------------------------------------------
      ps = p0

      !!!! modified by icon
      CALL init_velocity_field(deep, pertt, X, lon, lat, z, t, u, v)
      ! !-----------------------------------------------------
      ! !   Initialize velocity field
      ! !-----------------------------------------------------
      ! inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0)
      ! bigU = g / aref * K * inttau2 * inttermU * t
      ! if (deep .eq. 0) then
      !   rcoslat = aref * cos(lat)
      ! else
      !   rcoslat = (z + aref) * cos(lat)
      ! end if
  
      ! omegarcoslat = omegaref * rcoslat
      
      ! u = - omegarcoslat + sqrt(omegarcoslat**2 + rcoslat * bigU)
      ! v = 0.d0
  
      ! !-----------------------------------------------------
      ! !   Add perturbation to the velocity field
      ! !-----------------------------------------------------
  
      ! ! Exponential type
      ! if (pertt .eq. 0) then
      !   u = u + evaluate_exponential(lon, lat, z)
  
      ! ! Stream function type
      ! elseif (pertt .eq. 1) then
      !   u = u - 1.d0 / (2.d0 * dxepsilon) *                       &
      !       ( evaluate_streamfunction(lon, lat + dxepsilon, z)    &
      !       - evaluate_streamfunction(lon, lat - dxepsilon, z))
  
      !   v = v + 1.d0 / (2.d0 * dxepsilon * cos(lat)) *            &
      !       ( evaluate_streamfunction(lon + dxepsilon, lat, z)    &
      !       - evaluate_streamfunction(lon - dxepsilon, lat, z))
      ! end if
  
      !-----------------------------------------------------
      !   Initialize surface geopotential
      !-----------------------------------------------------
      phis = 0.d0
  
      !-----------------------------------------------------
      !   Initialize density
      !-----------------------------------------------------
      rho = p / (Rd * t)
  
      !-----------------------------------------------------
      !   Initialize specific humidity
      !-----------------------------------------------------
      if (moist .eq. 1) then
        eta = p/p0
  
        if (eta .gt. moisttr) then
          q = moistq0 * exp(- (lat/moistqlat)**4)          &
                      * exp(- ((eta-1.d0)*p0/moistqp)**2)
        else
          q = moistqs
        end if
  
        ! Convert virtual temperature to temperature
        t = t / (1.d0 + Mvap * q)
  
      else
        q = 0.d0
      end if
  
      !-----------------------------------------------------
      !   Initialize virtual potential temperature
      !-----------------------------------------------------
      thetav = t * (1.d0 + 0.61d0 * q) * (p0 / p)**(Rd / cp)
  
    END SUBROUTINE baroclinic_wave_test

  !-----------------------------------------------------------------------
  !  Init velocity field
  !-----------------------------------------------------------------------

    SUBROUTINE init_velocity_field(deep, pertt, X, lon, lat, z, tempv, u, v)
      INTEGER , INTENT(IN)  :: deep      ! Deep (1) or Shallow (0) test case
      INTEGER , INTENT(IN)  :: pertt     ! Perturbation type
      REAL(r8), INTENT(IN)  :: X         ! Earth scaling parameter
      REAL(r8), INTENT(IN)  :: lon       ! Longitude (radians)
      REAL(r8), INTENT(IN)  :: lat       ! Latitude (radians)
      REAL(r8), INTENT(IN)  :: z         ! Altitude (m)
      REAL(r8), INTENT(IN)  :: tempv     ! virtual temperature (K)
      REAL(r8), INTENT(OUT) :: u         ! zonal velocity
      REAL(r8), INTENT(OUT) :: v         ! meridional velocity
  
      ! local
      REAL(r8) :: aref, omegaref
      REAL(r8) :: inttermU
      REAL(r8) :: bigU
      REAL(r8) :: rcoslat
      REAL(r8) :: omegarcoslat
      REAL(r8) :: rratio
      REAL(r8) :: T0, constH, constC, scaledZ, inttau2
  
      !-----------------------------------------------------
      !   Initialize velocity field
      !-----------------------------------------------------
      aref = a / X
      omegaref = omega * X
  
      T0 = 0.5d0 * (T0E + T0P)
      constH = Rd * T0 / g
      constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
      scaledZ = z / (B * constH)
      inttau2 = constC * z * exp(- scaledZ**2)
  
      !--------------------------------------------
      !    radius ratio
      !--------------------------------------------
      if (deep .eq. 0) then
         rratio = 1.d0
      else
         rratio = (z + aref) / aref;
      end if
  
      inttermU = (rratio * cos(lat))**(K - 1.d0) - (rratio * cos(lat))**(K + 1.d0)
      bigU = g / aref * K * inttau2 * inttermU * tempv  !t
      if (deep .eq. 0) then
         rcoslat = aref * cos(lat)
      else
         rcoslat = (z + aref) * cos(lat)
      end if
  
      omegarcoslat = omegaref * rcoslat
  
      u = - omegarcoslat + sqrt(omegarcoslat**2 + rcoslat * bigU)
      v = 0.d0
  
      !-----------------------------------------------------
      !   Add perturbation to the velocity field
      !-----------------------------------------------------
  
      ! Exponential type
      if (pertt .eq. 0) then
         u = u + evaluate_exponential(lon, lat, z)
  
         ! Stream function type
      elseif (pertt .eq. 1) then
         u = u - 1.d0 / (2.d0 * dxepsilon) *                       &
              ( evaluate_streamfunction(lon, lat + dxepsilon, z)    &
              - evaluate_streamfunction(lon, lat - dxepsilon, z))
  
         v = v + 1.d0 / (2.d0 * dxepsilon * cos(lat)) *            &
              ( evaluate_streamfunction(lon + dxepsilon, lat, z)    &
              - evaluate_streamfunction(lon - dxepsilon, lat, z))
      end if
  
    END SUBROUTINE init_velocity_field
  
  !-----------------------------------------------------------------------
  !    Calculate pointwise pressure and temperature
  !-----------------------------------------------------------------------
    SUBROUTINE evaluate_pressure_temperature(deep, X, lon, lat, z, p, t)
  
      INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
  
      REAL(r8), INTENT(IN)  :: &
                  X,          & ! Earth scaling ratio
                  lon,        & ! Longitude (radians)
                  lat,        & ! Latitude (radians)
                  z             ! Altitude (m)
  
      REAL(r8), INTENT(OUT) :: &
                  p,          & ! Pressure (Pa)
                  t             ! Temperature (K)
  
      REAL(r8) :: aref, omegaref
      REAL(r8) :: T0, constA, constB, constC, constH, scaledZ
      REAL(r8) :: tau1, tau2, inttau1, inttau2
      REAL(r8) :: rratio, inttermT
  
      !--------------------------------------------
      ! Constants
      !--------------------------------------------
      aref = a / X
      omegaref = omega * X
  
      T0 = 0.5d0 * (T0E + T0P)
      constA = 1.d0 / lapse
      constB = (T0 - T0P) / (T0 * T0P)
      constC = 0.5d0 * (K + 2.d0) * (T0E - T0P) / (T0E * T0P)
      constH = Rd * T0 / g
  
      scaledZ = z / (B * constH)
  
      !--------------------------------------------
      !    tau values
      !--------------------------------------------
      tau1 = constA * lapse / T0 * exp(lapse * z / T0) &
           + constB * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
      tau2 = constC * (1.d0 - 2.d0 * scaledZ**2) * exp(- scaledZ**2)
  
      inttau1 = constA * (exp(lapse * z / T0) - 1.d0) &
              + constB * z * exp(- scaledZ**2)
      inttau2 = constC * z * exp(- scaledZ**2)
  
      !--------------------------------------------
      !    radius ratio
      !--------------------------------------------
      if (deep .eq. 0) then
        rratio = 1.d0
      else
        rratio = (z + aref) / aref;
      end if
  
      !--------------------------------------------
      !    interior term on temperature expression
      !--------------------------------------------
      inttermT = (rratio * cos(lat))**K &
               - K / (K + 2.d0) * (rratio * cos(lat))**(K + 2.d0)
  
      !--------------------------------------------
      !    temperature
      !--------------------------------------------
      t = 1.d0 / (rratio**2 * (tau1 - tau2 * inttermT))
  
      !--------------------------------------------
      !    hydrostatic pressure
      !--------------------------------------------
      p = p0 * exp(- g / Rd * (inttau1 - inttau2 * inttermT))
  
    END SUBROUTINE evaluate_pressure_temperature
  
  !-----------------------------------------------------------------------
  !    Calculate pointwise z and temperature given pressure
  !-----------------------------------------------------------------------
    SUBROUTINE evaluate_z_temperature(deep, X, lon, lat, p, z, t)
      
      INTEGER, INTENT(IN)  :: deep ! Deep (1) or Shallow (0) test case
  
      REAL(r8), INTENT(IN)  :: &
                  X,          & ! Earth scaling ratio
                  lon,        & ! Longitude (radians)
                  lat,        & ! Latitude (radians)
                  p             ! Pressure (Pa)
  
      REAL(r8), INTENT(OUT) :: &
                  z,          & ! Altitude (m)
                  t             ! Temperature (K)
  
      INTEGER :: ix
  
      REAL(r8) :: z0, z1, z2
      REAL(r8) :: p0, p1, p2
  
      z0 = 0.d0
      z1 = 10000.d0
  
      CALL evaluate_pressure_temperature(deep, X, lon, lat, z0, p0, t)
      CALL evaluate_pressure_temperature(deep, X, lon, lat, z1, p1, t)
  
      DO ix = 1, 100
        z2 = z1 - (p1 - p) * (z1 - z0) / (p1 - p0)
  
        CALL evaluate_pressure_temperature(deep, X, lon, lat, z2, p2, t)
        ! print*, ABS((p2 - p)/p)
        IF (ABS((p2 - p)/p) .lt. 1.0d-13) THEN
          EXIT
        END IF
  
        z0 = z1
        p0 = p1
  
        z1 = z2
        p1 = p2
      END DO
  
      z = z2
      ! print*, p, " pa => ",z," m" 
      CALL evaluate_pressure_temperature(deep, X, lon, lat, z, p0, t)
  
    END SUBROUTINE evaluate_z_temperature
  
  !-----------------------------------------------------------------------
  !    Exponential perturbation function
  !-----------------------------------------------------------------------
    REAL(r8) FUNCTION evaluate_exponential(lon, lat, z)
  
      REAL(r8), INTENT(IN)  :: &
                  lon,        & ! Longitude (radians)
                  lat,        & ! Latitude (radians)
                  z             ! Altitude (meters)
  
      REAL(r8) :: greatcircler, perttaper
  
      ! Great circle distance
      greatcircler = 1.d0 / pertexpr &
        * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
  
      ! Vertical tapering of stream function
      if (z < pertz) then
        perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
      else
        perttaper = 0.d0
      end if
  
      ! Zonal velocity perturbation
      if (greatcircler < 1.d0) then
        evaluate_exponential = pertup * perttaper * exp(- greatcircler**2)
      else
        evaluate_exponential = 0.d0
      end if
  
    END FUNCTION evaluate_exponential
  
  !-----------------------------------------------------------------------
  !    Stream function perturbation function
  !-----------------------------------------------------------------------
    REAL(r8) FUNCTION evaluate_streamfunction(lon, lat, z)
  
      REAL(r8), INTENT(IN)  :: &
                  lon,        & ! Longitude (radians)
                  lat,        & ! Latitude (radians)
                  z             ! Altitude (meters)
  
      REAL(r8) :: greatcircler, perttaper, cospert
  
      ! Great circle distance
      greatcircler = 1.d0 / pertr &
        * acos(sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon))
  
      ! Vertical tapering of stream function
      if (z < pertz) then
        perttaper = 1.d0 - 3.d0 * z**2 / pertz**2 + 2.d0 * z**3 / pertz**3
      else
        perttaper = 0.d0
      end if
  
      ! Horizontal tapering of stream function
      if (greatcircler .lt. 1.d0) then
        cospert = cos(0.5d0 * pi * greatcircler)
      else
        cospert = 0.d0
      end if
  
      evaluate_streamfunction = &
          (- pertu0 * pertr * perttaper * cospert**4)
  
    END FUNCTION evaluate_streamfunction
  
  END MODULE dcmip_baro_wave_test_mod